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We analyze two-dimensional (2D) random systems driven by a symmetric Levy stable noise which, 
under the sole influence of external (force) potentials $(1), asymptotically set down at Boltzmann- 
type thermal equilibria. Such behavior is excluded within standard ramifications of the Langevin 
approach to Levy flights. There, the action of a conservative force field ~ — V$(x) stands for an 
explicit reason for the emergence of an asymptotic invariant probability density function (pdf). 
However the latter cannot be represented in the Boltzmann form p*(x) ~ exp[— <E>(ie)] and the 
thermal equilibrium concept appears to be alien to Langevin-modeled Levy flights. In the present 
paper we address the response of Levy noise not to an external conservative force field, but directly 
to its potential That is explicitly encoded in non-symmetric jump transition rates of the 

master equation for the pdf p(x,t). We prescribe a priori the target pdf p, in the Boltzmann 
form ~ exp[— $(a;)] and next select the Levy noise of interest. Given suitable initial data, this 
allows to infer a reliable path-wise approximation to a true (albeit analytically beyond the reach) 
solution of the pertinent master equation, with the property p(x,t) — > p*(x) as time t goes to 
infinity. No explicit path-wise description has been so far devised for such thermally equilibrating 
random motion. To reconstruct random paths of the underlying stochastic process we resort to 
numerical methods, where long jumps of the Levy stable processes are statistically significant, but 
are truncated to become amenable to simulation procedures. We create a suitably modified version 
of the time honored Gillespie's algorithm, originally invented in the chemical kinetics context. A 
statistical analysis of generated sample trajectories allows us to infer a surrogate pdf dynamics which 
consistently sets down at a pre-defined target pdf. We pay special attention to the response of the 
2D Cauchy noise to an exemplary locally periodic "potential landscape" &(x),x € R 2 . 

I. INTRODUCTION 

Various random processes in real physical systems admit a simplified description based on stochastic differential 
equations. Then, there is a routine passage procedure from microscopic random variables to macroscopic (statistical 
ensemble, mean field) data, like e.g. the time evolution of an associated probability density function (pdf) which is 
a solution of a deterministic transport equation. A paradigm example is so-called Langevin modeling of diffusion- 
type and jump-type processes. The presumed microscopic model of random dynamics is provided by the Langevin 
(stochastic) equation, which additively decomposes into a (Newtonian by origin) drift and purely random (perturbing 
noise) term. Its direct consequence is the Fokker-Planck equation for an associated probability density function (pdf), 
and 0. 

As a necessary pre-requisite for our further discussion, let us discuss a transformation of the Fokker-Planck equa- 
tion into the Schrodinger-type (generalized diffusion) equation, often employed in the theoretical framework of the 
Brownian motion. Here, the Langevin equation, the induced Fokker-Planck equation and its Schrodingcr -type image 
are dynamically equivalent and describe the same diffusion- type process. This is not the case if one turns over to 
jump-type processes. 

For clarity of arguments, let us consider the Langevin equation for a one-dimensional diffusion process in an external 
conservative force field F(x) = —dV(x)/dx in the form % = F(x) + \/2vb{t), where b(t) stands for the normalized 
white noise: (b(t)) = 0, (b(t')b(i)) = 6(t — t') and the mass parameter is scaled away. The corresponding Fokker-Planck 
equation for the probability density function p(x, t) reads 

d t p = isAp-V(Fp) (1) 

and, in the confining regime, is known to enforce the existence of an asymptotic invariant pdf, p(x,t) — > p*(x) as 
t — > 00, in the explicit Boltzmann form exp(— 3>/2), where $(2;) = V{x)/v. 

By means of a standard substitution p(x,t) = rf;(x,t) exp[— V(x)/2u], the Fokker-Planck equation can be 
transformed into a generalized diffusion equation for an auxiliary function ?p(x,t). This Schrodinger-type equation 
(no imaginary unit i) reads 

d t i> = vAip - V(x)i> (2) 

where V{x) = \0^ ; + Vf) and F = F{x). 

By reintroducing a normalization constant (divide and multiply by a suitable number Z 1 / 2 in the factorization 
formula for p{x, t)), we can rewrite p(x, i) in the form p{x,t) = ^!{x, t)p 1 J 2 {x), where p^J 2 = Z^ 1 / 2 cxp(— < f ) /2) while 
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iff = Z 1 / 2 ^. Clearly, iff(x,t) — > pV 2 {x) as t goes to infinity. Moreover, we can rewrite the semigroup potential as 
follows: V{vApl /2 )/pl /2 . 

The transformation of (1) into (2) cannot be adopted to Levy jump-type processes, where the Langevin and 
Schrodinger-type (semigroup) modeling are known to be incompatible. Moreover, the Eliazar-Klafter no go statement, 
[H , disconnects the Langevin- modeled Fokker-Planck equation for any Levy-stable noise 

x = b{x) + A"(t) d t p = -V(6 • p) - A|A|^ /2 p . (3) 

from the very notion of the Boltzmann thermal equilibrium. 

However, the thermal equilibrium notion remains a valid concept within an immediate Levy transcript of the 
semigroup dynamics (2) (e.g. replace vA by — A|A| M / 2 ): 

d t V = -A|A|^/ 2 * - V*, (4) 

see e.g. [3, H|, where we assume that ^f(x,t) asymptotically sets down at a square root p 1 J 2 {x) of a well defined pdf 
p*. The semigroup potential V(x) follows from the compatibility condition: 

V = ~ A 172 • ( 5 ) 

p* 

In this particular context, while adopting a multiplicative decomposition of the time-dependent pdf p(x, t) — > p*(x): 

p{x,t) = V(x,t)pl /2 (x), (6) 

a novel fractional generalization of the Fokker-Planck equation governing the time evolution of p(x, t) has been 
introduced in Refs. see also [EH, 03, to handle systems that are randomized by symmetric Levy-stable drivers 

and may asymptotically set down at Boltzmann-typc equilibria under the influence of external potentials (thus not 
Newtonian forces anymore). 

The pertinent Fokker-Planck type equation, whose origin has been discussed before in a number of papers, (6j-[Io|. 
has the familiar master equation form, presently reproduced in the explicit 2D form: 

dtp(x, y) = jj [w^(x, y\u, v)p{u, v) - w<j,(u, v\x, y)p{x, y)]v^(x, y, dx, dy), (7) 

A 

Here, anticipating the effectiveness of numerical routines to be described in below, from the start we impose cut-offs 
upon the size of jumps to be accounted for during the simulations (large jump sizes remain statistically significant if 
view of the involved Levy distribution) 

A = {(x,y)eR 2 ; ef < \x - tt| < ef A e\ < \y - v\ ^ e%}, 
w$(x, y\u, v) = exp[($(u, v) - $(x, y))/2], 

and the Levy measure v^{x, y, dx, dy) in R 2 is given by 

2T((2 + p) 12) 1 _ 

^(x >y ,dx,dy)= ^ [r( _ M/2) | j^^y^m dxd y = M^y) dxd y^ ms(o,2). (9) 

It is the quantity w^, (x,y\u, v)v^ (x, y) which has an interpretation of the jump transition rate from the point (u, v) G R 2 
to another point {x,y) G R 2 . The potential function <f}(x,y) can be chose quite arbitrarily. However, we need to 
secure a L : (i? 2 ) normalization of the target pdf p*(x,y) ~ exp(— $>(x, y)). We note that exp(— $(x, y)) becomes a 
genuine stationary solution of Eq. ([7]) once we let e\' v —t and e^ v — > oo. 



II. GILLESPIE'S ALGORITHM: FINE TUNING IN 2D. 



The Gillespie's algorithm has been originally constructed to give a dynamical picture of a finite chain of chemical 
reactions, (Til Il2j |. There, switches between different chemical reaction channels can be re-interpreted as jumps 
between points in a finite state space. Since the number of allowed chemical channels is finite, a serious modification 
of an original algorithm must be created to account for very large (virtually infinite) state space we need to consider 
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FIG. 1: Gaussian target; statistics of 100000 random paths: surrogate pdf evolution and its OXY projection for a) t — 0.2, b) 
t — 0.7, c) t = 15. The subfigure d) refers to the asymptotic pdf (|14[) . All trajectories were started form the origin (0, 0). 



in connection with Levy nights. As an example, a jump process analog of chemical reaction channels could comprise 
(take R 1 for a while) all jumps form a fixed point xq to any point in the set [xq — £2, %o — ei] U [xq + ci,xq + e^]. Since 
we aim at numerical simulations, it is obvious that all admissible jumps cannot be computer-generated. Therefore, 
we must restrict considerations to a large number of suitably selected "representative" jumps of a truncated jump 
process ruled by the truncated Levy distribution of jumps. 

Basic tenets of the modified Gillespie's algorithm, fine tuned to account for Levy jumps in R 1 have been described 
in detail elsewhere, [llj]. Presently, we shall give a brief outline of the algorithm that is capable to account for Levy 
flights in R 2 . We mimic previously devised ID steps, (ljjj . while adopting them to the 2D situation. 



4 



Eq. (J7J can be rewritten as follows 

dtp{x,y) = // Cp[w$(x,y\£ + x,ri + y)p(£ + x ) T) + y)-w l f > (€ + x,ri + y\x,y)p(x,y)]d4dri, (10) 



A' 

where 

A' = {(£,77) G R 2 ; < |£| < e * A e\ ^ \ V \ < ef}, 

^ 2^sin( 7 r / z/2)[r(( / , + 2)/2)] 2 (11) 

P 7T 2 (£ 2 + t;2)( M +2)/2 

The algorithm outline reads as follows: 

(i) Set time t = and the point of origin for jumps, (xo,yo) € R ■ 

(ii) Create the set of all admissible jumps from (xo,yo) to (xo + £,yo + t?) , compatible with the transition rate 
w<l>(x,y\x + t;,y + ri). 

(iii) Evaluate 

Wi(x ,y )= J J C tl w tj> {£, + x ,'n + yo\xo,yo)d^dr], W 2 (x , yo) = J J C^w^ + x ,r] + yo\x ,yo)d£dr), 

-el -el -el e \ 



P P a 



(12) 

W 3 (x ,y ) = j j C,j,w$(£ + xo,r} + yo\xo,yo)d£dri, W 4 (x , y ) = j J C^w^ + x ,r] + y \x ,y )d£dr], 

and W^(a;o,yo) = Wi(x ,2/o) + W 2 (x ,yo) + W 3 {x ,y ) + Wi[x ,y ). 

(iv) Using a random number generator draw p x ,Py G [0, 1] from a uniform distribution 

(v) Employing above p x and p y and the identities 



/ / C^w^ + x , rj + y \x , y )d-q +J C lL + x , r) + y \x ,yo)dri \ d£ = p^W(x Q , y ), p x < Wi 2 (x , yo)/W(x , yo); 

-2 \-4 4 J 

k./-? 4 \ 

Wi2(x ,y ) + J J C^Mi^ff + x , V + yo\xo, yo)dr\ + / C^Wj,^ + x , r) + y \x , yo)dri I = p x W(x , yo), 1? W12 (%o, Vo)/W(x , y ); 

-t \-4 4 J 

/ / C M iu^(C + x , v + Vo\x , yo)d6 + J C^w^it + xq, 17 + yo\x , yo)d£ d)j = p y W(x , y ), p y < W 13 (x , yo)/W(x , yo); 



y / - e l E 2 \ 

W 13 (xo,y )+J J C^w^d + x , n + yo\x Q ,yo)d£, + j C m uj^(C + x , V + yo\x a , yo)d£ \ dr) = p y W(x ,y ), p y ^ W 13 (x , yo)/W(x a , yo), 

4 V'2 =f / 

(13) 

where Wi 2 (x 0l y ) = Wi(x ,y ) + W 2 (x ,y ) i Wi 3 (x , y ) = Wi(x ai y a ) + W 5 {xo,yo), find b x and b y corresponding 
to the "transition channel" (xo,yo) ~^ (bx,b y ) ■ 

(vi) Draw a new number q 6 (0, 1) from a uniform distribution. 

(vii) Reset time label s t = t + At where : At = — lnq/W(xo, yo). 

(viii) Reset (xo,yo) to a new location (xq = xq + b x ,yo = y + b y ), considered as a reference point for the iterative 
procedure. 

(ix) Return to step (ii) and repeat the procedure anew. 

Remark: All precautions respected in ID need to be respected in 2D as well (c.f. Comment 1 in [HI]). The 
following jump size bounds (integration boundaries) were adopted for exemplary numerical procedures to be 
described in below: sf — e\ = £\ = 0.001 i e 2 — e\ — e 2 — 1, provided one keeps in memory our convention not to 
reproduce the x, y index in e^'| anymore. 



III. STATISTICS OF RANDOM PATHS IN 2D: CASE STUDIES OF THE PDF EVOLUTION. 



Our main purpose in the present section is to analyze the response of free Levy noise to generic external potentials. 
That has been previously found to be encoded in jump transition rates. Now we address the same problem in a 
path-wise fashion. We shall faithfully follow the outlined random path generation procedure. Once suitable path 
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FIG. 2: y — projection of the previous pdf data (100000 trajectories) for running time instants a) t — 0.2. b) t = 0.7, c) 
t = 15, d) a projection of the target pdf is depicted. 
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a) b) 

FIG. 3: Time evolution of a) the first moment (average r = \/x 2 + y 2 value) , b) second moment (mean value of r 2 ). All 
100000 trajectories originate from (0,0). 



ensemble data are collected, we shall verify whether statistical (ensemble) features of generated random trajectories 
are compatible with those predicted by the master equation ([7]). This includes a control of an asymptotic behavior 
p(x, t) — > p*(x) with x <E R 2 , when t — > oo. 



A. Harmonic confinement: Gaussian target 



Let us prescribe an asymptotic invariant pdf p* to be in a 2D Gaussian form: 

P*{xiy) = -e-* 2 -y\ (14) 

7T 

As an exemplary source of random noise we assume the Cauchy driver, i.e. 2D Levy-stable noise with the stability 
index p = 1. Accordingly, the jump transition rate reads: 

I e -? 2 /2-x£-ri 2 ~yri 

C 1 w 4 ,(Z + x,r ] + y\x,y) = — ^ 2+t?3)3 /a • ( 15 ) 

To generate sample paths of the process, we need first to evaluate integrals of (|15[) over suitable integration volumes. 
If the integration volume comprises x,y which are close to the jump size boundaries ±£i, one develops a numerator 
of an expression into Taylor series and keeps terms up to the quadratic one. Errors induced by such approximation 
procedure are marginal. On the other hand, if x, y are far away from ±£i integrals (|15[) are amenable to standard 
evaluation methods (like e.g. Simposon's one). To be more explicit in this respect, let £12 = 0.05. If \a\, \b\, \c\, \d\ ^ £12, 
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then 



b d 



,-i 2 /2-x^- V 2 -y V 
(£ 2 + 7? 2 ) 3 / 2 



b d 



-d£ drj 



ie+v 2 ) 3/2 



2 



v 2 ) 



(16) 



where 



b d 



d£ drj 



ie + v 2 ) 3/2 



d$ = 



Vb 2 + d 2 Va 2 + d 2 Vb 2 + c 2 Va 2 + c 2 



b d 



x£d£ drj 

(e+v 2 ) 3/2 



= x In 



a c 
b d 



eVeTd? eVtFTc 2 

'(d + Vb 2 + d 2 )(c+Va 2 + c 2 ) 
{c + Vt) T T^){d+Va T Td 2 ) 



bd 



ad 



be 



ac 



a c 
b d 



y-qdid-q _ ^ I (b + Vb 2 + d 2 )(a + Va 2 + c 2 ) 
V n \(a + V^ 2 Td 2 )(b+ V^T^) 



ie+v 2 ) 3/2 

xy&ldS, df] 

(£2 + ^2)3/2 



xy 



(V« 2 + d 2 + \/b 2 +c 2 - ^b 2 +d 2 - Va 2 + c 2 ), 



a c 
b d 



x 2 - 1 £ 2 d£ d?7 

"2 (£2 + ^2)3/2 _ 2 



x 2 - 1 / „ b+^b 2 + d 2 
dm 



+ Va 2 + d 2 



a c 



y 2 - 1 ?/ 2 d£ dry y 2 - 1 / , , /d + Vb 2 + d 2 
ii6 = / — - — — ^tttt = — - — I bin 



2 (£2 + ^2)3/2 



c + y/b 2 + c 2 



— cln 



a In 



& + \Jb 2 + c 2 
a + Va 2 + c 2 



d + Va 2 + d 2 
c + \/a 2 + c 2 



(17) 



Numerical routines were written in terms of C-codes, 14 1. 

In Fig. 1 we have depicted the statistical data inferred from 100 000 trajectories, for three running time instants 
t = 0.2, t = 0.7, t = 15, together with the asymptotic expression (fT4|) . The right-hand-side column depicts projections 
of those data upon the the OXY plane. A substanitial increase of the analyzed enesemble data (like e.g. 300000, 
500000 or 10 6 ) is merely a matter of the simulation time span and adds nothing inspiring to the obtained behavior. 

It is clear that the surrogate pdf evolution consistently goes towards an invariant asymptotic pdf (compare subfigures 
c) and d)). We note a lowering and flattening of the maximum around (0, 0), in consistency with the ultimate target 
outcome, whose height is 1/ir as follows from p*(0,0) = l/ir. Visually accessible in homogeneities of circular shapes 
in OXY projections b) i c), are a consequence of still relatively low number (100 000) of sample paths data and 
approximations involved in evaluating involved integrals. The shape of OXY projection a) is a consequence of the 
initial data choice. Our problem has a radial symmetry. Therefore in Fig. 2 we depict a projection of the trajectory 
induced data upon the y ~ plane. The projection shows as well a consistent convergence towards the target pdf. 

An additional control method for the path-wise inferred pdf evolution, addresses the time evolution and an asymp- 
totic behavior of the pdf moments < R(t) > and < R 2 (t) >. Here < R(t) > is the mean distance of (x, y) points of a 
trajectory form the origin (0, 0) at the running time instant t, while < R 2 (t) > is a mean square distance from (0, 0). 
In view of 



<R> C 



1 



Vx^T 



0.886, 



< r2 >».= - H(x 2 + y 2 )e- x2 - 



(18) 



the < R(t) > dynamics should set down at \/tt/2, while this of < R 2 (t) > at 1. Fig. 3 depicts the evolution of 
< R(t) > and < R 2 (t) >, inferred from the simulated sample of 100000 jumping paths. 

The observed convergence < R(t) >^ 1/n and < R 2 {t) >— > 1 validates the number generator choice, we have used 
to arrive at sample jumping paths. 
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FIG. 4: Cauchy target; surrogate pdf evolution inferred from 100000 sample trajectories and the OXY projection for running 
times a) t = 0.2, b) t = 3, c) t = 500. A subfigure d) refers to the target Cauchy pdf I19p . In the course of simulations all 
trajectories were started from (0,0). 



B. Logarithmic confinement: 2D Cauchy target. 



We consider the target pdf p* in the 2D Cauchy form: 
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FIG. 5: The projection of the previous pdf data on the y — plane, at time instants a) t = 0.2, b) t = 3, c) t = 500, d) target 
pdf projection. 



Like previously, we take the Cauchy driver, fi = 1, as a reference Levy stable noise. Accordingly: 

„2 , „,2 \ 3/4 



Cimv(£ + x,rj + y\x,y) 



1 



1 



i + x z + y 2 



2n (e + v 2 ) 3/2 \ i + + 2 + {y + v) 2 

Proceeding like in the Gaussian case, for small £ and 77, Eq. (|20j) can be approximated by 

( - 3x 3y 



2tt (£ 2 + ?y 2 ) 3 / 2 



1 



2(1 + a: 2 +?y 2 ) * 2(l + a; 2 + y 2 ) 



(20) 



(21) 



where terms linear in £ i 77 were preserved. The result can be analytically integrated term after term by employing 
Eqs. (17).. 

Accumulated trajectory data have been analyzed to produce Fig. 4, where a surrogate pdf evolution is displayed. 
Fig. 5 reproduces the y = projection of the obtained pdf data. If compared with the ID case analyzed in Ref. [l3| 
a convergence to Boltzmannian equilibrium (target pdf) is substantially slowed down. 



C. Locally periodic confinement in R 



We consider target pdf whose Boltzmanian exponent is locally periodic (within a finite rectangle) and almost entirely 
localized within a finite spatial area due to harmonically confining tails of the potential: 



P*(x,y) = q 




\x\ ^ 2 i \y\ < 2 

\x\ > 2 i \y\ > 2 

W >2i \y\ ^2 

\x\ <2i \y\ >2, 



(22) 



The normalization constant C actually reads C = 13.5921. 

Subsequently adopted numerical integration routines heavily rely on the experience gained during our previous case 
studies. In Fig. 6 we report the surrogate pdf evolution, inferred from 100000 sample trajectories. A convergence 
rate to the asymptotic (target) pdf is satisfactory, although a reasonable agreement with the target data has been 
achieved for relatively large running time values, here t = 300. 

We are aware of the fact that the number of 100000 trajectories may be considered as too small and not sufficiently 
representative sample. Our tentative 300000 paths data do not show significant qualitative changes in the obtained 
evolution picture. 

We should mention that there are significant statistical fluctuations to be kept under control. They become are 
very conspicuous if the number of involved trajectories gets significantly lowered by imposing constraints (like e.g. 
various spatial projections). All 1 trajectory data, after being gathered, are safely stored in the computer memory. 
Therefore we can get access to any conceivable and more detailed statistical picture of what is going on, even if the 
outcome is hampered by significant random deviations from the reference (target) pdf data.. 

A sample of such fluctuating data is provided in Fig. 7, where we have considered projections of the surrogate pdf 
data upon planes y = 0, y = 1/4 and y = x at time t = 300. We have set them in a direct comparison with respective 
target pdf (|T^|) data. 



-4 -2 2 4 

X 

FIG. 6: Locally periodic target; surrogate pdf evolution and its OXY projection inferred from 100000 trajectories at running 
time instants a) t — 0.2, b) t = 3, c) t = 300. The subfigure d) depicts the asymptotic pdf ()22[) . All trajectories have been 
started from (0, 0). 



IV. OUTLOOK 

We have taken into consideration jump-type processes which can not be handled by standard stochastic differential 
equation methods (e.g. the Langevin modeling, where a conspicuous motion "tendency" quantified by an additive 
drift term, can be unambiguously isolated from the noise contribution) . Existing popular algorithms cannot provide a 
direct numerical simulation of sample paths of such non-standard processes. In the present paper, we have proposed 
a working method to generate stochastic trajectories (sample paths) of a random jump-type process that avoids any 
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FIG. 7: Projections of the surrogate pdf data at t — 300 upon planes a) y = 0, b) y = 1/4, c) y = x. An "asymptote" refers 
to the target pdf projections 



reference to a stochastic differential equation. An additional gain of that procedure is that we are in fact capable of 
reliably approximating the time evolution p(x, i) of a true (typically not available in a closed analytic form) solution 
of the master equation. 

To this end we have modified the Gillespie algorithm, [ll|, E3 > normally devised for sample paths generation if the 
transition rates refer to a finite number of states of a system. The essence of our modification is that we take into 
account the continuum of possible transition rates, thereby changing the finite sums in the original Gillespie algorithm 
into integrals. The corresponding procedures for stochastic trajectories generation have been changed accordingly. 

In other words, here we are able (i) to extract the background sample paths of a jump process and (ii) to infer a 
reliable approximation of an actual (analytically unavailable) solution of the master equation (7)- (8). We emphasize 
once more here, that we have focused on those jump-type processes that cannot be modeled by any stochastic 
differential equation of the Langevin type. 

Although heavy-tailed Levy stable drivers were involved in the present considerations, we have clearly confirmed 
that a large variety of stationary target distributions is dynamically accessible for each particular fi <= (0, 2) Levy driver 
choice. That variety comprises not only the standard Gaussian pdf, casually discussed in relation to the Brownian 
motion (e.g. the Wiener process), but the whole non-Gaussian family, associated with the Levy stable conceptual 
imagery . 

Among heavy-tailed distributions, we have paid attention to the Cauchy pdf which can stand for an asymptotic 
target for any (i € (0, 2) driver, provided a steering environment (e.g. "potential landscape") is properly devised. In 
turn, the Cauchy driver, while excited in a proper environment, may lead to an asymptotic pdf with an arbitrarily 
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large number of moments, the previously mentioned Gaussian case being included. . An example of the locally 
periodic environment has been considered as a toy model for more realistic physical systems. Our major hunch are 
strongly inhomogeneous "potential landscapes", modeled by relatively smooth potentials. We note that a radically 
extreme variant could be random potentials of Ref . Q ■ 

In connection with the master equation which was our departure point let us stress that, even if various mean 
field data are available in experimentally realizable systems, it is of vital interest to gain some knowledge about the 
microscopic dynamics (random paths) realized by the random system under consideration. The detailed analysis 
of sample path data with a focus on their specific features like e.g. ergodicity, mixing or lack of those properties, 
deserves a separate analysis. These goals can be achieved as well within the present simulation framework. It suffices 
to re-analyze the path-wise data we have collected and stored in the trajectory generation process. 
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